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Abstract 

In compressible turbulent combustion/nonequilibrium flows, the constructions of numerical schemes for 
(a) stable and accurate simulation of turbulence with strong shocks, and (b) obtaining correct propaga- 
tion speed of discontinuities for stiff reacting terms on “coarse grids” share one important ingredient 
- minimization of numerical dissipation while maintaining numerical stability. Here “coarse grids” 
means standard mesh density requirement for accurate simulation of typical non-reacting flows. This 
dual requirement to achieve both numerical stability and accuracy with zero or minimal use of numeri- 
cal dissipation is most often conflicting for existing schemes that were designed for non-reacting flows. 
The goal of this paper is to relate numerical dissipations that are inherited in a selected set of high order 
shock-capturing schemes with the onset of wrong propagation speed of discontinuities as a function of 
stiffness of the source term and the grid spacing. 


Motivation and Results 

To make accurate predictions in compressible turbulent combustion/nonequilibrium flows, one has 
to deal with the equations that describe time-dependent non-equilibrium effects, combustion, advanced 
thermodynamic models, and magnetic fields (MHD). Numerical simulation is challenging because of 
the conflicting requirements for numerical methods to be accurate enough to resolve the small scales of 
(magnetized) turbulence but robust enough to handle shock waves without generating spurious numeri- 
cal noise. In addition, the different physics models have different time scales that, when underresolved, 
might interact numerically to produce erroneous results. Furthermore, the appearance of the source 
terms in modeling non-equilibrium flow problems containing finite-rate chemistry or combustion poses 
additional numerical difficulties beyond that for solving non-reacting turbulent flows. The so-called stiff 
source term problem [5] is a well-known example. For stiff reactions it is well known that the wrong 
propagation speed of discontinuities occurs due to the under-resolved numerical solutions in both space 
and time. Schemes to improve the wrong propagation speed of discontinuities for systems of stiff react- 
ing flows remain a challenge for algorithm development [7] . 

In addition to the minimization of numerical dissipation while maintaining numerical stability in 
compressible turbulence with strong shocks, Yee & Sjogreen and Yee & Sweby [9, 10, 8] discussed a 
general framework for the design of such schemes. Yee & Sjogreen [12], Wang et al. [7] and references 
cited therein present their recent progress on the subject. In [13], a short overview of this recent progress 
is given. Two very important numerical challenges are “Stiffness and Nonlinearity of Source Terms”. 

The objective of the present paper is to a gain a deeper understanding on the behavior and perfor- 
mance of four high order shock-capturing schemes for two representative stiff detonation wave prob- 
lems. The test cases consist of the Arrhenius ID Chapman- Jouguet (C-J) detonation wave [4, 6] and a 
2D Heaviside detonation wave [1]. These are the same two test cases considered in [7]. The consid- 
ered four schemes are the fifth-order WENO, “WEN05”, the [7] newly developed subcell resolution 
version of WEN05, “WEN05/SR”, the Yee et al. nonlinear filter version of WEN05, “WEN05fi” 
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Figure 1: Density (mass fraction of the unburnt gas; 50 grid points) and pressure (800 grid points) comparison in the vicinity 
of the discontinuity among four high order shock-capturing methods for the C-J detonation problem, Arrhenius case at t = 1.8 
using 50 and 800 uniform grid points. 


and the Durcros split version of WEN05fi, “WEN05fi+split” [11, 12]. For the temporal discretization, 
the classical fourth-order Runge-Kutta method (RK4) is used. Although WEN05/SR, WEN05fi and 
WEN05fi+split use the same dissipative portion of WEN05 as the basic numerical dissipation model, 
the mechanism to further limit the inherited amount of WEN05 dissipation from one grid point to the 
next are different among the schemes. See the aforementioned references for the development of these 
schemes. 

WEN05/SR [7] is a newly developed modified fractional step method which solves the convection 
step and reaction step separately. In the convection step any high order shock-capturing method can be 
used. In the reaction step an ODE solver is applied, but with the computed flow variables in the shock 
region modified by the Harten subcell resolution idea [3]. 

WEN05fi is the filter version of WEN05. On the first stage a full time step by RK4 is performed. 
For this stage the sixth-order central spatial base scheme is used. On the second stage the solution is fil- 
tered by the dissipative portion of WEN05 in conjunction with a wavelet flow sensor. The wavelet flow 
sensor indicates the locations where shock-capturing dissipations are needed and leaves the remaining 
region free of numerical dissipation contamination. WEN05fi+split is WEN05fi applied to the Ducros 
et al. split form of the governing equation [2] before the application of WEN05fi. The Ducros et al. split 
form is a preprocessing step to condition the governing equation(s) before the application of high order 
central schemes. This preprocessing step improves numerical stability and is widely used in numerical 
modeling and simulation of turbulent flows. 

The comparison of the performance of the four schemes is largely based on the degree that each 
method captures the correct location and jump size of the stiff reaction front for coarse grids. It is 
remarked that, in order to resolve the sharp stiff reaction zone, sufficiently many grid points in this zone 
are still needed. The behavior of these schemes in the vicinity of a sharp reaction zone with several 
levels of grid refinement will be briefly touched upon. Figures 1 and 2 show the comparison among 
WEN05, WEN05/SR, WEN05fi and WEN05fi+split for the the Arrhenius ID Chapman- Jouguet (C- 
J) detonation wave [4, 6] and a 2D Heaviside detonation wave [1]. The reference solution for the ID stiff 
detonation wave case is computed by the regular WEN05 scheme with 10, 000 uniform grid points. The 
initial values consist of totally burnt gas on the left-hand side and totally unburnt gas on the right-hand 
side. The reference solution for the 2D stiff detonation wave case is computed by standard WEN05 with 
2000 x 400 grid points. The chemical reaction for this 2D problem is modeled by the Heaviside form. 
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Figure 2: ID cross-section of density at t = 1.7 x 10 7 by four high order shock-capturing methods for the 2D detonation 
problem using 200 x 40 uniform grid points with left figure in the vicinity of the discontinuity. 


The grid density for all computed solutions is considered to be very coarse for the subject combustion 
test cases. From the result, WEN05/SR and WEN05fi+split are able to obtain the correct shock speed 
with similar accuracy whereas this is not the case for WEN05 WEN05fi using the same coarse grids. 
Using its original form [12] without further modification, the accuracy of WEN05fi+split is nearly as 
good as the newly developed high-order finite difference schemes with subcell resolution. The results 
indicate the influence of numerical dissipation on the behavior of the computed solution around the 
stiff reaction front. In the final paper, we will examine the subcell resolution version of WEN05fi and 
WEN05fi+split and their seventh-order counterparts. 
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